www.gusucode.com > Matlab在化学工程中的应用 > Matlab在化学工程中的应用/实用化工计算机模拟－Matlab在化学工程中的应用/Examples/Chapter 5/PDEs2DS_CrankNicolson.m
function PDEs2DS_CrankNicolson
% 用有限差分求解固定床反应器二维拟均相稳态模型（二维稳态PDE方程组）
%
% Author: HUANG Huajiang
% Copyright 2003 UNILAB Research Center,
% East China University of Science and Technology, Shanghai, PRC
% $Revision: 1.0 $ $Date: 2002/05/02 $
global n F G rc dz M M1 b1 b2 C2 C3 C4
% Increment of r and z
n = 5+1; % node number in r direction is 5
m = 100;
dr = 0.01; % m
dz = 0.0585; % m
% Parameters
Ramda = 0.45; % W/(m K)
G = 2500; % kg/(m^2 h)
Cp = 2.18; % kJ/(kg K)
dH = 140e3; % J/mol
rho = 1440; % kg/m^3
a = 0.05; % radius of reactor, m
u0c0 = 0.069/(pi*a^2); % u0*c0 obtained from u0*c0*A=0.069 kmol/h
Cp1 = 1.0e3; % 烟道气比热，J/(kg K)
F = 130/3600; % kg/s
% Equation coefficient（方程的系数）
a1 = Ramda/(G*Cp*1000)*3600; % m
b1 = dH*rho/G/Cp; % Note of unit accordance：dH*rho/G/Cp rc = [K/m]，rc=[kmol/(kg h)]
a2 = 0.000427; % a2 = D2/mu
b2 = rho/u0c0;
C1 = 2*pi*a*Ramda/(F*Cp1); % formula (12)
C2 = 4*dr/(dz*C1); % formula (33)
M = a1*dz /dr^2; % formula (17)
M1 = a2*dz /dr^2; % formula (19)
C3 = M*C2/2*(1+0.5/(n-1))-(1+M);
C4 = M*C2/2*(1+0.5/(n-1))-(1-M);
% Initializing variables, at the inlet of the reactor
% (在反应器入口处即z=0或j=1时初值)
T(1:n-1) = 873;
T(n) = 893;
x(1:n) = 0;
X0 = [T;x];
Tresult(1,:) = X0(1,:);
xresult(1,:) = X0(2,:);
% Calculate F(i,j),G(i,j) and rc(i,j) at z=0 (initial value)
% 计算对应于j=1时的F(i,1)和G(i,1)及反应速度rc
FGrc(T,x);
for j=2:m
X = fsolve(@TxEquations,X0); % Solve nonlinear equations
T = X(1,:);
x = X(2,:);
Tresult(j,:) = T
xresult(j,:) = x
xa(j,:) = sum(2*xresult(j,2)+4*xresult(j,3)+6*xresult(j,4)...
+8*xresult(j,5)+5*xresult(j,6))/25
if xa(j)>0.45 % When the average conversion > 0.45, stop calculating and exit the loop
break
end
% Calculate F(i,j),G(i,j) and rc(i,j) for the next iteration
% 计算F(i,j)和G(i,j)及反应速度rc(i,j)供下次迭代使用
FGrc(T,x)
end
z = dz.*[0:j-1];
r = dr.*[0:n-1];
% When the conversion rate of ethylbenzene is 45%, the reactor length required is:
z = z'
L = spline(xa,z,0.45)
% Plot the results
surf(r,z,Tresult) % 反应管轴径向温度分布
xlabel('r (m)')
ylabel('z (m)')
zlabel('T (K)')
figure
plot(z,xa) % 平均转化率沿管长的分布图
xlabel('z (m)')
ylabel('x_a_v')
figure
surf(r,z,xresult) % 轴径向平均转化率分布
xlabel('r (m)')
ylabel('z (m)')
zlabel('x')
% --------------------------------------------------------------------------
function f = ReactionRate(T,x)
% Calculate the reaction rate（计算反应速度）
k = 0.027*exp(0.021*(T-773));
f = 15100*exp(-11000./T).*((1-x)./(11+x)-1.2*x.^2./k./(11+x).^2);
% --------------------------------------------------------------------------
function f = FGrc(T,x)
% Calculate F(i), G(i)
global F G rc n dz M M1 b1 b2 C2 C3 C4
rc = ReactionRate(T,x); % Calculate the reaction rate
F(1) = ((1-2*M)*T(1)+2*M*T(2)-b1*dz/2*rc(1))/(1+2*M); % formula (28)
G(1) = ((1-2*M1)*x(1)+2*M1*x(2)+b2*dz/2*rc(1))/(1+2*M1); % formula (29)
i = (2:5);
var1 = 1-0.5./(i-1);
var2 = 1+0.5./(i-1);
F(i) = ( M/2*( var1.*T(i-1)+var2.*T(i+1) ) + (1-M)*T(i) -b1*dz/2*rc(i) )/(M+1); % formula (22)
G(i) =( M1/2*( var1.*x(i-1)+var2.*x(i+1) ) + (1-M1)*x(i)+b2*dz/2*rc(i) )/(M1+1); % formula (23)
F(n) = ( -M*T(n-1) +C4*T(n)+b1*dz/2*rc(n) )/ C3; % formula (35)
G(n) = ( M1*x(n-1)+(1-M1)*x(n)+b2*dz/2*rc(n) )/(1+M1); % formula (31)
% --------------------------------------------------------------------------
function f = TxEquations(X)
% Define the linear equation system composed of (20)、(21) 、(26)、(27)、(30) and (34).
global n F G rc dz M M1 b1 b2 C2 C3
T = X(1,:)
x = X(2,:)
fT(1) = F(1) + (2*M*T(2)-b1/2*dz*rc(1) )/(1+2*M) - T(1); % formula (26)
fx(1) = G(1) + (2*M1*x(2)+b2/2*dz*rc(1) )/(1+2*M1) - x(1); % formula (27)
for i=2:n-1
var1(i) = (1-0.5/(i-1));
var2(i) = (1+0.5/(i-1));
fT(i) = F(i)+ ( M/2*(var1(i)*T(i-1)+var2(i)*T(i+1))-b1/2*dz*rc(i) )/(M+1) - T(i); % formula (20)
fx(i) = G(i)+ ( M1/2*(var1(i)*x(i-1)+var2(i)*x(i+1))+b2/2*dz*rc(i) )/(M1+1) - x(i); % formula (21)
end
fT(n) = F(n) + ( -M*T(n-1)+b1/2*dz*rc(n) )/C3 -T(n); % formula (34)
fx(n) = G(n) + ( M1*x(n-1)+b2/2*dz*rc(n) )/(1+M1) - x(n); % formula (30)
f = [fT;fx];